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Due to the stringent design requirement for aerospace or nuclear 
structural components, considerable research interests have been gener- 
ated on the development of constitutive models for representing the 
inelastic behavior of metals at elevated temperatures. In particular, 
a class of unified theories (or viscoplastic constitutive models) have 
been proposed [1-10] to simulate material responses such as cyclic 
plasticity, rate sensitivity, creep deformations, strain hardening or 
softening, etc. This approach differs from the conventional creep and 
plasticity theory in that both the creep and plastic deformations are 
treated as unified time-dependent quantities. Although most of visco- 
plastic models give better material behavior representation, the associ- 
ated constitutive differential equations have stiff regimes which present 
numerical difficulties in time-dependent analysis. In this connection, 
appropriate solution algorithm must be developed for viscoplastic 
analysis via finite element method. 

In the past, inelastic finite element structural analyses were per- 
formed largely based on the classical concept of creep and plasticity 
[11-14], Recently, some attempts have been made to incorporate a speci- 
fic type of viscoplastic theories- into finite element codes [15-20] for 
structural analysis. In this paper, three integration schemes are 
implemented into a nonlinear finite element program [21] to study their 
numerical efficiency pertaining to finite element analysis. Moreover, 
four viscoplastic models, namely, those due to Walker, Miller, Krieg- 
Swearingen-Rohde, and Robinson, were implemented into a finite element 
program for nonlinear analysis. A general implementation procedure is 
outlined in the paper. 
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VISCOPLASTIC THEORIES 


The basic assumption embodied in viscoplastic theories is the 
unified treatment of inelastic strain; i.e. no distinction is given to 
creep and plastic deformations. In addition, both elastic and inelastic 
strains are considered to be present at all stages of loading and un- 
loading processes. The unique feature of such treatment, as compared to 
classical theories, is that the yield condition Is not explicitly in- 
volved. Consequently, the computational algorithm for complex loading 
history can be much simplified. In the context of small deformation, 
viscoplastic models may be written in the following general form 

* C * * 

a 3 D e - (2GeI + s (a + 3G)yT) (1) 

( 2 ) 

(3) 

(4) 

T 3 Temperature; 

0 3 Stress vector; 

K 3 Drag stress; 

X,G = Lame contant; 
el 3 Inelastic strain vector; 

Y 3 Linear thermal expansion coefficient; 
f 3 Inelastic strain rate function; 
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Time devative; 

Total strain vector; 

Back stress vector; 
Elasticity matrix; 
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he,,^ = Hardening functions for back and drag stresses, respectively; 
•"a^k 3 Recovery functions for back and drag stresses, respectively. 

Eqs. (1-4) represent a complete set of viscoplastic constitutive 
equations wherein the following assumptions are invoked in the extention 
from uniaxial case to the three-dimensional case, namely, i) isotropic 
material, ii) incompressible inelastic strain, and iii) linear bulk be- 
havior. Eq. (1) defines stress rate to be proportional to elastic strain 
rate while Eq. (2) states the functional dependence of inelastic strain 
rate on applied stress, temperature, and state variables. Furthermore, 
Eqs. (3-4), so-called evolutional equations, are generally constructed in 
hardening/recovery form such that the net effect of two antagonistic 
mechanisms uniquely determines the growth rate of state variables a and 
k. 

Although the mathematical expressions of viscoplastic models pro- 
posed by various researchers differ in their detailed descriptions, they 
do however portray several common phenomena: i) Initial linear elastic 
behavior wherein the inelastic effect is negligible and then nonlinear 
response afterwards, ii) strain-rate sensitivity, iii) time-dependent 
creep and relaxation, iv) cyclic hardening or softening, v) creep recov- 
ery, vi ) creep plasticity interactions, and vii) Bauschinger's effect. 

NUMERICAL INTEGRATION SCHEMES 

For finite element applications, it is useful to choose an appro- 
priate integration scheme for handling the nonlinear viscoplastic 
equations concerned. Krieg (22) pointed out the exi stance of numerical 
stiff regions in viscoplastic formulation together with a discussion of 
potential difficulties. The stiffness of the equations originates from 
the nonlinear relationship assumed in Eq. (1) and the hardening/recovery 
form in evolutional equations. Formal definition of the stiffness of a 
set of differential equations can be found in Ref. (23) where the measure 
of "stiffness" is given in terms of the spectra of eigenvalues obtained 
from the Jacobian matrix of associated equation system. 

Numerical approaches intended for integrating stiff differential 
equations have been developed by a number of researchers. Among them 
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Gear's method is the most famous one. Although Gear's package has been 
used quite effectively in solving one-dimensional constitutive equations, 
(1) it is not suitable for large scale finite element analysis simply 
because its solution procedure is of a multistep nature. When employed 
in finite element analysis, this method usually requires a large amount 
of storage in order to follow the deformation history of the material. 
For this reason, one-step method is much preferable. 

For the purpose of discussion, the constitutive equations are re- 
written as follows 


y * f (y t) ( 5 ) 

rw /v /v * 

where y represents stress, inelastic strain and state variables while f 
denotes nonlinear functions. One-step method for solving inelastic rate 
problems in the field of finite element has been investigated by several 
researchers (24-26). In a broad sense, it can be written in terms of 
one-parameter (0) family of implicit algorithm (the - method) as follows. 


y n +l = Yn 4t n [ (1-0) fn + 9 ^n+1 3 ^ 

~ ~ ~ ~ 

where At n = t n +i - t n is the n-th time step size and 0 is an integration 
parameter which has the range of (0,1). In Eq. (6) it is assumed that a 
numerical solution at the beginning of time step n is known, the solution 

at the end of the step is to be sought. 

The simplest integration scheme is the explicit forward Euler scheme 
corresponding to 0 = 0. It is an explicit scheme since the solution at 
time t n+ i is completely determined from conditions existing at time t n . 
Therefore, in the forward Euler method, the solution at time t n+ i is 
approximated by 

Yn+1 = yn + At n f n 

/V ^ ^ 

When this method is employed in solving stiff equations, very small step 
size must be used in order to obtain stable and accurate solutions. 

On the other hand, the case 0 - 1/2 results in the so-called impli- 
cit trapezoidal scheme which is also widely known as Crank-Nicholson rule 
in the context of linear differential equations. Then 
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(8) 


*n+l s *n + 

r* r* 


At n 

— (£ n + £ n +l) 


Note that f n+ j = f(y n +i, t n+ i) Is unknown. Nonlinear implicit equation 

is best solved by the Newton-Raphson iteration. To this end, Eq. (8) is 
rewritten in the form 

1 1 1 

F n+1 a yn+1 ' yn " ^n( f n+l + F n )/2 a 0 (9) 


The right superscript "i " denotes iteration number, 
known, Newton-Raphson iteration gives 





Since y n and f n are 


( 10 ) 


Rearranging Eq. (10) yields 
-i 

,Xl 
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Defining 


1+1 (i+1) i 

Ay = y - y 

~n+l ~n+l ~n+l 


(ID 


( 12 ) 


and performing differentiation, one obtains 


At 3f ' . 

r — _!!Ii i+i i “ i 

U - 2 1 J^n-Kl - Jt„ - i n+ l * 2 (f n * f„ +1 ) 

dy n+l 


(13) 


where the initial value of y may be obtained by an explicit scheme. 

n+1 


Eq. (13) stands for a linear system of equations for implicit trape- 
zoidal method. The system is readily solved by Gaussian elimination and 
backward substitution. If this method is employed in an analysis, the 
immediate question is: how one can determine whether the solution has 
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converged or not? 
for this purpose, 
y such that 


e 


In fact, several convergence criteria could be used 
One convenient way is to check the iterative value of 


ail! 2 

y 


<Tol 


(14) 


where || || = Euclidean norm 

Tol a A tolerance ratio 

Presently, the above criterion is employed to determine the convergence 
of a solution. 

Comparing Eqs. (7) & (13), it is apparent that the implicit trape- 
zoidal method requires not only much more functional evaluations but also 
solving a system of linear equations. As an alternative, the implicitness 
of f n+ i in Eq. (8) may be removed by using Talor series expansion, namely. 


fn+1 * ^n + Jn^yn+l 
~ • ~ * ~ 

(15) 

where 

J n = 3f n /3yn 

(16) 

Thus, Eq. (13) becomes 



[I - J n At/2]Ay n 3 At f n (17) 

22 22 ~ 

The above equation is referred as the explicit trapezoidal scheme since 
the solution is completely determined from the initial conditions. 

At this point, it is instructive to make some qualitative comparisons 
among the aforementioned numerical schemes. Comparing explicit trapezoi- 
dal scheme with forward Euler scheme reveals that they differ only in the 
expression J n At/2, i.e. the product of Jacobian matrix and half of step 
size. The addition of such matrix necessitates the solution be obtained 
by solving a system of simultaneous equations. Like implicit trapezoidal 
scheme, it also requires the evaluations of Jacobian matrix. Apparently, 
by including the extra term, the numerical behavior of the constitutive 
equations have become stabilized. In this context, J n At/2 essentially 
plays the role of a correcting factor. On the other hand, since no itera- 
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tlon Is Involved in the explicit trapezoidal scheme, it can be viewed as a 
starter of Implicit trapezoidal scheme. 

We consider another extreme case, i.e. 0 ■ 1, which is called impli- 
cit Euler scheme, 

y n +l 3 yn + at n f n+ i (18) 

** *** 

Hughes et al (24) demonstrated that for viscoplastic finite element anal; - 
sis one-parameter family of implicit algorithm is unconditionally stable 
when e > 1/2 while only conditionally stable otherwise. In recent years, 
various numerical schemes have been applied to viscoplastic problems 
(15-20,24-26). Some of the authors have discarded the explicit Euler 
method due to its numerical instability. However, the validity of this 
conclusion needs to be further explored. In Ref. (27), present investi- 
gators evaluated three numerical techniques for integrating the visco- 
plastic constitutive equations for a uniaxial state of stress. The 
schemes evaluated were: i) forward Euler method, ii) explicit trapezoidal 
method, and 111) implicit trapezoidal method with Newton-Raphson itera- 
tion method. Although implicit trapezoidal method with iteration appears 
to be more stable and accurate than the other methods even when the step 
size is considerable large, its suitability for finite element analysis 
must be re-assessed. 

In principle, inelastic analysis using finite element method con- 
sists of a sequence of Incremental process. Two most widely used 
approaches are the Initial strain and tangent stiffness methods. In con- 
sideration of the formulation presented in Eqs. (1-4), one finds that the 
initial strain method is the most natural way to handle viscoplastic 
models. The reason behind this will be elaborated below. 

In Eq. (1), we invoke an assumption that the strain increment is 
decomposed into elastic and inelastic components. Then, the inelastic 
part, which is governed entirely by Eqs. (1-4) at constitutive level, is 
converted into an equivalent load in the finite element formulation. 

Thus, we have 

kEau = (Ap 0 ) + (Ap e ) (19) 

M os* ** 

where 

K e s Elastic structural stiffness matrix, which may vary with 
the temperature 
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Au = Incremental nodal displacement vector 
Ap 0 * Incremental vector of applied load 

APe “ Incremental vector due to inelastic and thermal strains. 

In addition to the incremental procedure used for solving the global 
stiffness equations, a subincrementing technique Is employed to calculate 
the constitutive material matrix. That is, let At be the time increment 
for solving the global stiffness equations. Then At is sub-divided into 
smaller increments with a constant step size, at = At/m. Moreover, the 
number of subincrements can be determined by an automatic stepping pro- 
cedure for which an error measure is compared with a specified tolerance. 
Further discussion of this is given in [27]. 

COMPUTER IMPLEMENTATION 

With the constitutive relations and numerical integration schemes 
outlined in the previous sections, the next step Is to implement these 
relationships into a typical (general purpose) finite element program for 
Intended analysis. For this purpose, the related computer subroutines 
are written in the form of an Independent material module so that it can 
be easily interfaced with a finite element code. 

The calculation steps for a viscoplastic model can be summarized as 

follows: 

Step 1. Preset the strains, stresses, back stresses, inelastic strains, 
nodal temperatures, etc. transferred from the main program. 

Step 2. For non-isothermal condition, interpolate temperature at Gauss 
points from nodal temperatures. 

Step 3. Compute strain rate and temperature rate, and select step size 
of subincrements. 

Step 4. Interpolate temperature dependent material constants based on 
the average temperature at the mid-point of a time step. 

Step 5. Solve for the state variables from the constitutive equations 
using a subincrementing method with a selected integration 
technique. 

Step 6. Check for solution convergence and determine whether cut-back of 
step size is necessary. 
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Step 7. Update the stresses, strains, inelastic strains and other state 
variables, then return to the main program. 

EXAMPLE 

To demonstrate the utility of the finite element procedures, Robin- 
son's unified theory was applied to the analysis of a pressurized thick- 
walled cylinder which is restrained in its axial direction. Finite 
element mesh, its dimension and boundary condition are shown in Fig.l. 

The loading history consists of a 0.0028 hour ramp up to an internal 
pressure of 3.65 ksi followed by a hold period at that temperature for 
200 hours. Explicit trapezoidal method was employed for this example. 

Figs. 2 to 3 show the predicted redistribution of hoop and axial 
stress at several selected time following rapid pressurization, wherein 
zero time denotes the end of the loading ramp. As can be seen, while the 
internal pressure Is held constant, these stresses undergo variation 
exhibiting rapid redistribution followed by a steady-steady response. 

The tendency of approaching to a saturated state is apparent. 

Figs. 4 and 5 show the creep displacement at the outside wall of 
thick-walled cylinder using both Euler and explicit trapezoidal scheme 
with different time step sizes as well as number of subincrements. Solid 
line indicates the base-line solution using the explicit trapezoidal 
scheme. From the numerical experiments, it was found that the forward 
Euler Integration scheme with an automatic stepping and error control Is 
far more efficient In computation as compared to the explicit and impli- 
cit trapezoidal schemes. 


REFERENCES 

1. Miller, A.K., "An Inelastic Constitutive Model for Monotonic, 

Cyclic, and Creep Deformation" ASME Journal of Engineering Materials 
and Technology, Vol. 98, 1976, pp. 97-113. 

2. Kreig, R.D., Swearingen, J.C., and Rohde, R.W., "A Physically-based 
Internal Variable Model for Rate-Dependent Plasticity", Inelastic 
Behavior of Pressure Vessel and Piping Components, T.Y. Chang and 
E. Krempl Eds., ASME PVP-PB-028, 1978, pp. 15-28. 


195 



3. Walker, K.P., "Representation of Hastelloy-X Behavior at Elevated 
Temperature with a Functional Theory of Viscoplaticity", Paper 
presented at ASME/PVP Century 2 Energy Technology Conference, San 
Francisco, CA., Aug. 1980. 

4. Hart, E.W. , "Constitutive Relations for the Nonelastic Deformation 
of Metals", ASME Journal of Engineering Material and Technology, 

Vol . 98, 1976, pp. 193-202. 

5. Robinson, D.N., "A Unified Creep-Plasticity Model for Structural 
Metals at High Temperature" ORNL/TM-5969, Oak Ridge National Labora- 
tory, Oct. 1978. 

6. Cernocky, E.P., and Krempl , E. , "A Theory of Theromo-Viscoplasticity 
Based on Infinitesimal Total Strain", International Jouranl of Solid 
and Structures, Vol. 16, 1980, pp. 723-741. 

7. Cernocky, E.P., and Krempl, E., "A Nonlinear Uniaxial Integral Con- 
stitutive Equation Incorporating Rate Effects, Creep, and Relaxation" 
International Journal of Nonlinear Mechanics, Vol. 14, 1979, pp. 
183-389. 

8. Bodner, S.R., and PArtom, Y., "Constitutive Equations for Elastic- 
Viscoplastic Strain-Hardening Materials", ASME Journal of Applied 
Mechanics, Vol. 42, 1975, pp. 385-389. 

9. Chaboche, J.L., and Rousselier, G., "On the Plastic and Viscoplastic 
Constitutive Equations", Inelastic Analysis and Life Prediction in 
Elevated Temperature Design, Edited by G. Baylac, ASME PVP-Vol. 59, 
1982, pp. 33-35. 

10. Larsson, B., and Storakers, B., "A State Variable Interpretation of 
Some Rate-Dependent Inelastic Properties of Steel", ASME Journal of 
Engineering Materials and Technology, Vol. 100, 1978, pp. 395-401. 


196 



11. Snyder, M.D., and Bathe, K.J., "A Solution Procedure for Thermo- 
Elastic-Plastic and Creep Problems", Journal of Nuclear Engineering 
and Oesign, Vol . 64, 1981, pp. 49-80. 

12. Haisler, W.E., and Sander, O.R., "Elastic-Plastic-Creep-Large Strain 
Analysis at Elevated Temperature by the Finite Element Method", 
Computer and Structures, Vol. 10, 1979, pp. 375-382. 

13. Cyr, N.A., and Teter, R.D., "Finite Element Elastic-Plastic-Creep 
Analysis of Two-Dimensional Continuum with Temperature Dependence 
Material Properties", Computer and Structures, Vol. 3, 1973, pp. 
849-863. 

14. Levy, A., and Pifko, A.B., "On Computational Strategies for Problems 
Involving Plasticity and Creep", International Journal for Numerical 
Methods in Engineering, Vol. 17, 1981, pp. 747-771. 

15. Morjaria, M., and Mukherjee, S., "Finite Element Analysis of Time- 
Dependent Inelastic Deformation in the Presence of Transient Thermal 
Stresses", International Journal of Numerical Methods in Engineering 
Vol. 17, 1981, pp. 909-921. 

16. Robinson, D.N., and Swinderman, R.W., "Unified Creep-Plasticity 
Constitutive Equations for 2-1/4 Cr - 1 Mo Steel at Elevated Temper- 
ature", ORNL/TM-8444, Oak Ridge National Laboratory, 1983. 

17. Zirin, R.M., and Krempl, E., "A Finite Element Integration Method 
for the Theory of Viscoplasticity Based on Infinitesimal Total 
Strain", ASME Journal of Pressure Technology, Vol. 104, May 1982, 
pp. 130-136. 

18. Newman, M., Zaphir, Z., and Bodner, S.R., "Finite Element Analysis 
for Time Dependent Inelastic Material Behavior", Computer and 
Structures, Vol. 6, 1976, pp. 157-162. 


197 



19. Walker, K.P., "Research and Development Program for Nonlinear 
Structure Modeling with Advanced Time-Dependent Constitutive Rela- 
tionships", Report to NASA Lewis Research Center, NAS3-22055. 

20. Chang, T.Y., NFAP - A Nonlinear Finite Element Analysis Program, 

Vol . 1 and 2, Dept, of Civil Engineering, The University of Akron, 
Akron, Ohio, Oct. 1980. 

21. Krieg, R.D., "Numerical Integration of Some New Unified Plasticity- 
Creep Formulation", Transactions of the 4-th SMIRT Conference, San 
Francisco, CA., Vol. M, 1977, Paper M 6/4. 

22. van der Houwen, P.J., Construction of Integration Formular for 
Initial Value Problem, North-Hoi land Publishing Company, Amsterdam, 

1977. 

23. Hughes, T.J.R., "Unconditionally Stable Algorithms for Quasi-Static 
Elasto/Visco-plastic Analysis", Computer and Structures, Vol. 8, 

1978, pp. 169-173. 

24. Wiliam, K.J., "Numerical Solution of Inelastic Rate Processes", 
Computer and Structures, Vol. 8, 1978, pp. 511-531. 


25. Cormeau, I., " Numerical Stability in Quasi-Static Elasto/Visco- 

Plasticity", International Journal for Numerical Methods in Engine 
ering, Vol. 9, 1975, pp. 109-127. 



. 075 " 

Fig. 1 - A Thick Wal'l Cylinder Under Internal Pressure 
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Fig. 3 - Axial Stress Distribution 
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Fig. 4 - Creep Displacement at Outside Wall of the Cylinder 

by Various Integration Steps. 



Fig. 5 - Creep Displacement Predicted by Euler Method with Automatic Stepping 
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